Planck distribution (planck)

Planck distribution (planck)#

SciPy’s planck distribution is a discrete exponential distribution on the non‑negative integers.

[ \mathbb{P}(X=k)=(1-e^{-\lambda})e^{-\lambda k},\qquad k=0,1,2,\dots\quad (\lambda>0) ]

It is mathematically the same as a geometric distribution (counting failures before the first success) under the reparameterization (p = 1-e^{-\lambda}).

This is not the continuous “Planck’s law” distribution of photon wavelengths/energies; it is the discrete distribution of occupation numbers (counts) for a single mode.

Learning goals#

  • Classify the distribution and state its support/parameter space.

  • Derive the PMF/CDF and connect them to geometric / truncated Boltzmann.

  • Compute mean/variance/skewness/kurtosis, MGF/CF, and entropy.

  • Derive the likelihood and the MLE for (\lambda).

  • Implement NumPy-only sampling and validate it with Monte Carlo.

  • Use scipy.stats.planck for PMF/CDF/RVS and fitting via scipy.stats.fit.

Prerequisites#

  • Comfort with geometric series and basic calculus.

  • Familiarity with log, exp, and numerical-stability helpers like expm1 and log1p.

Notebook roadmap#

  1. Title & Classification

  2. Intuition & Motivation

  3. Formal Definition

  4. Moments & Properties

  5. Parameter Interpretation

  6. Derivations

  7. Sampling & Simulation

  8. Visualization

  9. SciPy Integration

  10. Statistical Use Cases

  11. Pitfalls

  12. Summary

import numpy as np

import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio

from scipy import stats

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=6, suppress=True)

SEED = 7
rng = np.random.default_rng(SEED)

import sys
import scipy
import plotly

print("Python:", sys.version.split()[0])
print("NumPy:", np.__version__)
print("SciPy:", scipy.__version__)
print("Plotly:", plotly.__version__)
print("Seed:", SEED)

Python: 3.12.9
NumPy: 1.26.2
SciPy: 1.15.0
Plotly: 6.5.2
Seed: 7

1) Title & Classification#

  • Name: planck (Planck discrete exponential distribution)

  • Type: Discrete

  • Support (base distribution): (k \in {0,1,2,\dots})

  • Parameter space (SciPy shape): (\lambda>0)

SciPy also allows a location shift loc (an integer shift for discrete distributions). With loc:

[ X = \texttt{loc} + K,\qquad K\in{0,1,2,\dots}. ]

We will mostly work with the base distribution (K) (i.e. loc=0).

Notation and convenient reparameterizations:

  • (K \sim \mathrm{Planck}(\lambda))

  • Define (q = e^{-\lambda} \in (0,1)) and (p = 1-q = 1-e^{-\lambda}\in(0,1))

Then the PMF is (\mathbb{P}(K=k)=p,q^k), i.e. a geometric distribution on ({0,1,2,\dots}).

2) Intuition & Motivation#

What this distribution models#

The Planck distribution models a count with an exponentially decaying tail:

  • large counts are possible, but their probability shrinks like (e^{-\lambda k})

  • the distribution is right-skewed (often strongly)

A key interpretation comes from physics:

  • Consider a single quantum harmonic oscillator mode in thermal equilibrium.

  • Each additional excitation (“photon”, “quantum”) costs energy proportional to the count.

  • The probability of seeing (k) excitations is proportional to (e^{-\lambda k}), and normalization yields the Planck PMF.

Typical real-world use cases#

  • Occupation numbers (counts) of bosonic modes in statistical physics.

  • Any count process where the probability of additional units decays exponentially.

  • A simple discrete alternative to exponential/gamma when data are integer-valued (e.g. binned waiting times).

Relations to other distributions#

  • Geometric distribution: exactly the same distribution under (p = 1-e^{-\lambda}).

  • Boltzmann distribution (scipy.stats.boltzmann): a truncated discrete exponential on ({0,\dots,N-1}). As (N\to\infty), it approaches the Planck distribution.

  • Exponential distribution: if (E\sim\mathrm{Exp}(\text{rate}=\lambda)), then (\lfloor E\rfloor\sim\mathrm{Planck}(\lambda)).

  • Negative binomial: sums of independent geometric/Planck variables produce negative binomial counts.

A useful qualitative property: memorylessness. For (q=e^{-\lambda}), [ \mathbb{P}(K\ge k+n\mid K\ge k)=q^n. ] Among discrete distributions on ({0,1,2,\dots}), this is characteristic of the geometric/Planck family.

3) Formal Definition#

Let (K \sim \mathrm{Planck}(\lambda)) with (\lambda>0). Define (q=e^{-\lambda}\in(0,1)) and (p=1-q).

PMF#

[ \mathbb{P}(K=k) = (1-e^{-\lambda})e^{-\lambda k} = p,q^k,\qquad k\in{0,1,2,\dots} ] and (\mathbb{P}(K=k)=0) otherwise.

A numerically stable way to compute (1-e^{-\lambda}) for small (\lambda) is (1-e^{-\lambda} = -\mathrm{expm1}(-\lambda)).

CDF#

Because this is a discrete distribution, the CDF is step-like. For real (x): [ F(x)=\mathbb{P}(K\le x)= \begin{cases} 0, & x<0\ 1 - e^{-\lambda(\lfloor x\rfloor+1)}, & x\ge 0. \end{cases} ]

Equivalently, for an integer (k\ge 0): [ F(k) = 1 - q^{k+1}. ]

Survival function#

For an integer (k\ge 0): [ \mathbb{P}(K>k) = q^{k+1} = e^{-\lambda(k+1)}. ]

def validate_lambda(lambda_: float) -> float:
    lambda_ = float(lambda_)
    if not (np.isfinite(lambda_) and lambda_ > 0.0):
        raise ValueError(f"lambda_ must be finite and > 0, got {lambda_!r}")
    return lambda_


def _is_integer_array(x: np.ndarray) -> np.ndarray:
    """True for integer dtype and for floats that are exact integers."""
    return np.equal(x, np.floor(x))


def planck_logpmf(k, lambda_: float):
    """Log-PMF of the Planck distribution at integer k >= 0."""
    lambda_ = validate_lambda(lambda_)
    k = np.asarray(k)

    out = np.full_like(k, -np.inf, dtype=float)
    mask = (k >= 0) & _is_integer_array(k)

    # log(1 - exp(-lambda_)) computed stably
    log_p = np.log(-np.expm1(-lambda_))
    out[mask] = log_p - lambda_ * k[mask]
    return out


def planck_pmf(k, lambda_: float):
    """PMF of the Planck distribution at integer k >= 0."""
    return np.exp(planck_logpmf(k, lambda_))


def planck_cdf(x, lambda_: float):
    """CDF of the Planck distribution for real x."""
    lambda_ = validate_lambda(lambda_)
    x = np.asarray(x, dtype=float)

    out = np.zeros_like(x, dtype=float)
    mask = x >= 0
    k = np.floor(x[mask])

    # F(x) = 1 - exp(-lambda_ * (floor(x) + 1))
    out[mask] = -np.expm1(-lambda_ * (k + 1.0))
    return out

# Quick sanity checks vs SciPy

lambda_ = 1.3
ks = np.arange(0, 15)

pmf_np = planck_pmf(ks, lambda_)
pmf_sp = stats.planck.pmf(ks, lambda_)

print("max |pmf numpy - pmf scipy|:", float(np.max(np.abs(pmf_np - pmf_sp))))

K = 200
cdf_tail = planck_cdf(K, lambda_)
print(f"P(K <= {K}) = {cdf_tail:.12f} (so remaining tail mass ~ {1-cdf_tail:.2e})")

max |pmf numpy - pmf scipy|: 6.938893903907228e-18
P(K <= 200) = 1.000000000000 (so remaining tail mass ~ 0.00e+00)

4) Moments & Properties#

Let (q=e^{-\lambda}\in(0,1)) and (p=1-q).

Mean, variance, skewness, kurtosis#

For (K\sim\mathrm{Planck}(\lambda)):

Quantity

Closed form

Mean

(\mathbb{E}[K]=\dfrac{q}{p}=\dfrac{1}{e^{\lambda}-1})

Variance

(\mathrm{Var}(K)=\dfrac{q}{p^2}=\dfrac{e^{\lambda}}{(e^{\lambda}-1)^2}=\mu(1+\mu))

Skewness

(\gamma_1 = \dfrac{1+q}{\sqrt{q}} = 2\cosh(\lambda/2))

Excess kurtosis

(\gamma_2 = q + q^{-1} + 4 = 2\cosh(\lambda) + 4)

(Full kurtosis is (3+\gamma_2). SciPy’s stats(..., moments='k') returns excess kurtosis.)

MGF / characteristic function#

The probability generating function (PGF) is [ G(z)=\mathbb{E}[z^K] = \frac{p}{1-qz},\qquad |qz|<1. ]

  • MGF (M(t)=\mathbb{E}[e^{tK}]) exists when (t<\lambda) and equals [ M(t)=\frac{p}{1-qe^{t}} = \frac{1-e^{-\lambda}}{1-e^{t-\lambda}},. ]

  • Characteristic function (\varphi(\omega)=\mathbb{E}[e^{i\omega K}]) exists for all real (\omega): [ \varphi(\omega)=\frac{p}{1-qe^{i\omega}}. ]

Entropy#

The (Shannon) entropy is [ H(K) = -\sum_{k\ge 0} \mathbb{P}(K=k)\log \mathbb{P}(K=k) = -\log(1-e^{-\lambda}) + \frac{\lambda}{e^{\lambda}-1}. ]

Memorylessness#

[ \mathbb{P}(K\ge k+n\mid K\ge k)=q^n\quad\text{for integers }k,n\ge 0. ]

def planck_mean(lambda_: float) -> float:
    lambda_ = validate_lambda(lambda_)
    return float(1.0 / np.expm1(lambda_))


def planck_variance(lambda_: float) -> float:
    mu = planck_mean(lambda_)
    return float(mu * (1.0 + mu))


def planck_skew(lambda_: float) -> float:
    lambda_ = validate_lambda(lambda_)
    q = np.exp(-lambda_)
    return float((1.0 + q) / np.sqrt(q))


def planck_excess_kurtosis(lambda_: float) -> float:
    lambda_ = validate_lambda(lambda_)
    q = np.exp(-lambda_)
    return float(q + 1.0 / q + 4.0)


def planck_entropy(lambda_: float) -> float:
    lambda_ = validate_lambda(lambda_)
    mu = planck_mean(lambda_)
    return float(-np.log(-np.expm1(-lambda_)) + mu * lambda_)


def planck_mgf(t, lambda_: float):
    """MGF M(t) = E[e^{tK}] for t < lambda_."""
    lambda_ = validate_lambda(lambda_)
    t = np.asarray(t, dtype=float)
    if np.any(t >= lambda_):
        raise ValueError("MGF exists only for t < lambda_")

    q = np.exp(-lambda_)
    p = -np.expm1(-lambda_)
    return p / (1.0 - q * np.exp(t))


def planck_cf(w, lambda_: float):
    """Characteristic function phi(w) = E[e^{i w K}] for real w."""
    lambda_ = validate_lambda(lambda_)
    w = np.asarray(w, dtype=float)
    q = np.exp(-lambda_)
    p = -np.expm1(-lambda_)
    return p / (1.0 - q * np.exp(1j * w))


lambda_ = 1.3

mu = planck_mean(lambda_)
var = planck_variance(lambda_)
skew = planck_skew(lambda_)
exkurt = planck_excess_kurtosis(lambda_)
ent = planck_entropy(lambda_)

mean_sp, var_sp, skew_sp, exkurt_sp = stats.planck.stats(lambda_, moments='mvsk')
ent_sp = stats.planck.entropy(lambda_)

print("mean:", mu, "| scipy:", float(mean_sp))
print("var:", var, "| scipy:", float(var_sp))
print("skew:", skew, "| scipy:", float(skew_sp))
print("excess kurt:", exkurt, "| scipy:", float(exkurt_sp))
print("entropy:", ent, "| scipy:", float(ent_sp))

# MGF spot-check (finite sum vs closed form)
t = 0.5
K_max = 2000
ks = np.arange(0, K_max + 1)
mgf_sum = np.sum(np.exp(t * ks) * planck_pmf(ks, lambda_))
mgf_closed = float(planck_mgf(t, lambda_))
print(f"MGF(t={t}) sum ~ {mgf_sum:.8f} | closed {mgf_closed:.8f}")

mean: 0.3746305205153175 | scipy: 0.3746305205153175
var: 0.5149785474168952 | scipy: 0.5149785474168951
skew: 2.437586605774912 | scipy: 2.4375866057749125
excess kurt: 7.941828460653257 | scipy: 7.941828460653257
entropy: 0.8052046593263178 | scipy: 0.8052046593263178
MGF(t=0.5) sum ~ nan | closed 1.32105769

/tmp/ipykernel_1040229/2480825993.py:71: RuntimeWarning:

overflow encountered in exp

/tmp/ipykernel_1040229/2480825993.py:71: RuntimeWarning:

invalid value encountered in multiply

5) Parameter Interpretation#

(\lambda) controls how quickly the tail decays.

  • Larger (\lambda) (\Rightarrow) faster decay (\Rightarrow) most mass near 0.

  • Smaller (\lambda) (\Rightarrow) heavier tail (mean and variance increase like (\approx 1/\lambda) and (\approx 1/\lambda^2) for small (\lambda)).

Useful equivalences:

  • (q=e^{-\lambda}) is the tail ratio: (\mathbb{P}(K\ge k) = q^k).

  • (p=1-e^{-\lambda}) is the geometric success probability.

So you can interpret (\lambda) as a discrete decay rate (analogous to the rate of an exponential distribution).

# Shape changes: PMF curves for different lambda_

lambdas = [0.2, 0.5, 1.0, 2.0]
ks = np.arange(0, 35)

fig = go.Figure()
for lam in lambdas:
    fig.add_trace(
        go.Scatter(
            x=ks,
            y=planck_pmf(ks, lam),
            mode="lines+markers",
            name=f"lambda_={lam}",
        )
    )

fig.update_layout(
    title="Planck PMF for different lambda_",
    xaxis_title="k",
    yaxis_title="P(K = k)",
)
fig.show()

# Mean / variance as a function of lambda_
lam_grid = np.linspace(0.05, 3.0, 200)
mu_grid = np.array([planck_mean(l) for l in lam_grid])
var_grid = np.array([planck_variance(l) for l in lam_grid])

fig2 = go.Figure(
    data=[
        go.Scatter(x=lam_grid, y=mu_grid, mode="lines", name="mean"),
        go.Scatter(x=lam_grid, y=var_grid, mode="lines", name="variance"),
    ],
    layout=go.Layout(
        title="Mean and variance vs lambda_",
        xaxis_title="lambda_",
        yaxis_title="value",
        legend_title="moment",
    ),
)
fig2.show()

6) Derivations#

Throughout, let (q=e^{-\lambda}\in(0,1)) and (p=1-q).

Expectation#

Start from the PMF (\mathbb{P}(K=k)=p q^k): [ \mathbb{E}[K] = \sum_{k=0}^{\infty} k,p q^k = p\sum_{k=0}^{\infty} k q^k. ] Using the standard geometric-series identity (\sum_{k\ge 0} k q^k = \frac{q}{(1-q)^2}) (for (|q|<1)): [ \mathbb{E}[K] = p\frac{q}{(1-q)^2} = \frac{q}{1-q}. ] Substituting (q=e^{-\lambda}) gives (\mathbb{E}[K] = \frac{1}{e^{\lambda}-1}).

Variance#

Compute (\mathbb{E}[K^2]) using (\sum_{k\ge 0} k^2 q^k = \frac{q(1+q)}{(1-q)^3}): [ \mathbb{E}[K^2] = p\frac{q(1+q)}{(1-q)^3} = \frac{q(1+q)}{(1-q)^2}. ] Then [ \mathrm{Var}(K)=\mathbb{E}[K^2]-\mathbb{E}[K]^2 = \frac{q}{(1-q)^2}. ]

Likelihood and MLE#

For observations (k_1,\dots,k_n\in{0,1,2,\dots}), the log-likelihood is [ \ell(\lambda) = \sum_{i=1}^n \log\big((1-e^{-\lambda})e^{-\lambda k_i}\big) = n\log(1-e^{-\lambda}) - \lambda\sum_{i=1}^n k_i. ] Differentiate and set to zero: [ \ell’(\lambda) = \frac{n}{e^{\lambda}-1} - \sum_{i=1}^n k_i = 0 \quad\Longrightarrow\quad \widehat{\lambda} = \log\Big(1+\frac{n}{\sum_i k_i}\Big)=\log\Big(1+\frac{1}{\bar{k}}\Big). ]

If (\bar{k}=0) (all observations are zero), the likelihood increases as (\lambda\to\infty), reflecting a near-degenerate distribution at 0.

def validate_sample(k):
    k = np.asarray(k)
    if k.ndim != 1:
        raise ValueError("sample must be 1D")
    if np.any(~np.isfinite(k)):
        raise ValueError("sample must be finite")
    if np.any(k < 0) or np.any(~_is_integer_array(k)):
        raise ValueError("sample must contain integers >= 0 (use loc-shifted model if needed)")
    return k.astype(int)


def planck_loglik(lambda_: float, sample) -> float:
    lambda_ = validate_lambda(lambda_)
    k = validate_sample(sample)

    n = k.size
    s = int(np.sum(k))

    # log L = n log(1 - exp(-lambda_)) - lambda_ * sum k_i
    return float(n * np.log(-np.expm1(-lambda_)) - lambda_ * s)


def planck_mle_lambda(sample) -> float:
    k = validate_sample(sample)
    m = float(np.mean(k))
    if m == 0.0:
        return float("inf")
    return float(np.log1p(1.0 / m))


# Example: recover lambda_ from simulated data
true_lambda = 1.1
sample = stats.planck.rvs(true_lambda, size=2000, random_state=rng)

lam_hat_closed = planck_mle_lambda(sample)
print("true lambda_:", true_lambda)
print("closed-form MLE:", lam_hat_closed)

# Compare to generic optimizer-based fit
fit_res = stats.fit(stats.planck, sample, bounds={"lambda": (1e-6, 10.0)})
print("scipy.stats.fit params:", fit_res.params)

true lambda_: 1.1
closed-form MLE: 1.1066768485048402
scipy.stats.fit params: FitParams(lambda_=1.1066767722935067, loc=0.0)

7) Sampling & Simulation#

NumPy-only sampler (inverse transform via tail)#

The survival function is (\mathbb{P}(K\ge k)=q^k=e^{-\lambda k}).

If (U\sim\mathrm{Unif}(0,1)), define [ K = \Big\lfloor \frac{-\log U}{\lambda} \Big\rfloor. ] Then for integer (k\ge 0): [ \mathbb{P}(K\ge k) = \mathbb{P}\left(\frac{-\log U}{\lambda}\ge k\right) = \mathbb{P}(U\le e^{-\lambda k}) = e^{-\lambda k}, ] so (K) has the desired Planck PMF.

This is also a nice conceptual bridge:

  • (E=-\log U/\lambda) is an (\mathrm{Exp}(\text{rate}=\lambda)) random variable

  • and (K=\lfloor E\rfloor)

def planck_rvs_numpy(lambda_: float, size: int, *, rng: np.random.Generator) -> np.ndarray:
    lambda_ = validate_lambda(lambda_)
    size = int(size)

    u = rng.random(size)
    # rng.random can (rarely) return exactly 0, which would map to +inf.
    u = np.maximum(u, np.nextafter(0.0, 1.0))

    return np.floor(-np.log(u) / lambda_).astype(int)


# Monte Carlo check
lambda_ = 0.9
n = 200_000
s = planck_rvs_numpy(lambda_, n, rng=rng)

print("MC mean:", float(np.mean(s)), "| theory:", planck_mean(lambda_))
print("MC var :", float(np.var(s)), "| theory:", planck_variance(lambda_))

MC mean: 0.68208 | theory: 0.6851177504050078
MC var : 1.1418368736 | theory: 1.1545040823250263

8) Visualization#

We’ll visualize:

  • the PMF (bars)

  • the CDF (step curve)

  • a Monte Carlo histogram compared to the theoretical PMF

lambda_ = 0.8
ks = np.arange(0, 30)

pmf = planck_pmf(ks, lambda_)
cdf = planck_cdf(ks, lambda_)

fig_pmf = go.Figure(
    data=[go.Bar(x=ks, y=pmf, name="PMF")],
    layout=go.Layout(
        title=f"Planck PMF (lambda_={lambda_})",
        xaxis_title="k",
        yaxis_title="P(K = k)",
    ),
)
fig_pmf.show()

fig_cdf = go.Figure(
    data=[go.Scatter(x=ks, y=cdf, mode="lines+markers", name="CDF")],
    layout=go.Layout(
        title=f"Planck CDF (lambda_={lambda_})",
        xaxis_title="k",
        yaxis_title="P(K ≤ k)",
    ),
)
fig_cdf.show()

# Monte Carlo vs theory
n = 80_000
samples = planck_rvs_numpy(lambda_, n, rng=rng)

# Empirical frequencies for 0..K_max
K_max = 25
counts = np.bincount(samples[samples <= K_max], minlength=K_max + 1)
emp_pmf = counts / n

fig_mc = go.Figure()
fig_mc.add_trace(go.Bar(x=np.arange(K_max + 1), y=emp_pmf, name="Empirical", opacity=0.6))
fig_mc.add_trace(
    go.Scatter(
        x=np.arange(K_max + 1),
        y=planck_pmf(np.arange(K_max + 1), lambda_),
        mode="lines+markers",
        name="Theory",
    )
)

fig_mc.update_layout(
    title=f"Monte Carlo PMF vs theory (n={n:,}, lambda_={lambda_})",
    xaxis_title="k",
    yaxis_title="Probability",
)
fig_mc.show()

9) SciPy Integration#

SciPy provides scipy.stats.planck as an rv_discrete distribution.

Common methods:

  • pmf, logpmf

  • cdf, sf

  • rvs

  • stats (mean/var/skew/excess kurtosis)

  • entropy

Fitting#

Many rv_discrete distributions (including planck) do not implement a .fit(...) method. In modern SciPy, you can fit discrete or continuous distributions with scipy.stats.fit.

For planck, the MLE for lambda_ also has a closed form (derived above), which is often preferable.

lambda_ = 1.0
rv = stats.planck(lambda_)  # frozen distribution (loc=0)

ks = np.arange(0, 6)
print("pmf:", rv.pmf(ks))
print("cdf:", rv.cdf(ks))
print("rvs:", rv.rvs(size=10, random_state=rng))

mean_sp, var_sp, skew_sp, exkurt_sp = rv.stats(moments='mvsk')
print("mean/var/skew/excess kurt:", float(mean_sp), float(var_sp), float(skew_sp), float(exkurt_sp))
print("entropy:", float(rv.entropy()))

# Fit lambda_ from data using scipy.stats.fit (fix loc=0 implicitly)
data = rv.rvs(size=2000, random_state=rng)
fit_res = stats.fit(stats.planck, data, bounds={"lambda": (1e-6, 10.0)})
print("fit params:", fit_res.params)

lam_hat_closed = planck_mle_lambda(data)
print("closed-form MLE:", lam_hat_closed)

pmf: [0.632121 0.232544 0.085548 0.031471 0.011578 0.004259]
cdf: [0.632121 0.864665 0.950213 0.981684 0.993262 0.997521]
rvs: [0 0 0 0 1 2 1 0 0 1]
mean/var/skew/excess kurt: 0.5819767068693265 0.9206735942077924 2.2552519304127614 7.086161269630487
entropy: 1.0406518522564083
fit params: FitParams(lambda_=0.9855410710546102, loc=0.0)
closed-form MLE: 0.9855410569229435

10) Statistical Use Cases#

Hypothesis testing (likelihood-ratio test)#

A simple test for (H_0: \lambda=\lambda_0) uses the likelihood ratio statistic [ \Lambda = 2\big(\ell(\widehat{\lambda}) - \ell(\lambda_0)\big). ] Under regularity conditions and large (n), (\Lambda) is approximately (\chi^2_1) under (H_0).

Bayesian modeling (Beta–Geometric, then transform)#

Using the geometric parameterization (\mathbb{P}(K=k)=p(1-p)^k), with a Beta prior (p\sim\mathrm{Beta}(\alpha,\beta)), the posterior is conjugate: [ p\mid k_{1:n} \sim \mathrm{Beta}\Big(\alpha+n,; \beta+\sum_i k_i\Big). ] Then (\lambda) is a deterministic transform: (\lambda = -\log(1-p)).

Generative modeling#

Because Planck is geometric:

  • it is a natural count likelihood for exponentially-tailed discrete data

  • sums and mixtures connect it to common families (e.g. negative binomial)

  • it can serve as a simple component in hierarchical count models

def planck_lrt(sample, lambda0: float):
    """Likelihood ratio test for H0: lambda_ = lambda0."""
    lambda0 = validate_lambda(lambda0)
    sample = validate_sample(sample)

    lam_hat = planck_mle_lambda(sample)
    ll_hat = planck_loglik(lam_hat, sample) if np.isfinite(lam_hat) else 0.0
    ll_0 = planck_loglik(lambda0, sample)

    stat = 2.0 * (ll_hat - ll_0)
    p_value = float(stats.chi2.sf(stat, df=1))
    return stat, p_value, lam_hat


# Hypothesis test example
true_lambda = 1.2
sample = stats.planck.rvs(true_lambda, size=1500, random_state=rng)

lambda0 = 0.9
stat, p_value, lam_hat = planck_lrt(sample, lambda0=lambda0)

print("true lambda_:", true_lambda)
print("H0 lambda_:", lambda0)
print("MLE lambda_:", lam_hat)
print("LRT stat:", stat)
print("p-value :", p_value)


# Bayesian example: Beta prior on p = 1 - exp(-lambda_)
alpha0, beta0 = 2.0, 2.0
k = validate_sample(sample)

n = k.size
s = int(np.sum(k))

alpha_post = alpha0 + n
beta_post = beta0 + s

n_draws = 20_000
p_draws = rng.beta(alpha_post, beta_post, size=n_draws)

# Transform to lambda = -log(1 - p)
lambda_draws = -np.log1p(-p_draws)

fig = px.histogram(
    lambda_draws,
    nbins=60,
    title=f"Posterior over lambda_ (Beta prior on p, alpha={alpha0}, beta={beta0})",
    labels={"value": "lambda_"},
)
fig.add_vline(x=true_lambda, line_dash="dash", line_color="black", annotation_text="true")
fig.add_vline(x=lam_hat, line_dash="dash", line_color="red", annotation_text="MLE")
fig.show()

true lambda_: 1.2
H0 lambda_: 0.9
MLE lambda_: 1.1951784870591609
LRT stat: 99.78071204285334
p-value : 1.702401395301426e-23

11) Pitfalls#

  • Invalid parameters: (\lambda\le 0) is not valid.

  • Numerical cancellation for small (\lambda): computing (1-e^{-\lambda}) as 1 - np.exp(-lambda_) loses precision when (\lambda) is small.

    • Prefer -np.expm1(-lambda_) and np.log(-np.expm1(-lambda_)).

  • Huge means for small (\lambda): (\mathbb{E}[K]=1/(e^{\lambda}-1)\approx 1/\lambda) blows up as (\lambda\to 0).

    • Monte Carlo estimates may need very large sample sizes to stabilize.

  • Degeneracy as (\lambda\to\infty): the distribution collapses toward 0; higher standardized moments become ill-conditioned.

  • Fitting with location shifts: if your data are supported on ({c,c+1,\dots}), you need a loc=c shift (or subtract (c)) before using the base formulas.

12) Summary#

  • planck is a discrete exponential distribution on ({0,1,2,\dots}) with PMF (\mathbb{P}(K=k)=(1-e^{-\lambda})e^{-\lambda k}).

  • It is exactly a geometric distribution under (p=1-e^{-\lambda}), and has the memoryless tail (\mathbb{P}(K\ge k)=e^{-\lambda k}).

  • Closed forms are available for moments, MGF/CF, entropy, and the MLE (\widehat{\lambda}=\log(1+1/\bar{k})).

  • A simple NumPy-only sampler is (K=\lfloor -\log(U)/\lambda\rfloor) with (U\sim\mathrm{Unif}(0,1)).

  • In SciPy, use scipy.stats.planck for evaluation/sampling and scipy.stats.fit (or the closed-form MLE) for fitting.

Contents